#ifndef _SYSTEM_WS_H_
#define _SYSTEM_WS_H_
#include<string>
#include"base.h"
#include"orbital.h"
#include"system_ho.h"
#include"help.h"
#include"constants.h"
using std::string;
class System_Tz_WS:public System_Tz<WS_Orbital_Tz,complexd>
{
 public:
  typedef System_Tz_HO System_TzType;
  typedef WS_Orbital_Tz OrbitalType;
  typedef complexd DataType;
  typedef TwoBodyMatElement<DataType> TwoBodyMatEle;
  const static int MatEledim;//< dimension of TwoBodyMatEle.

  System_Tz_WS(){}
 System_Tz_WS(const string & _OneBodyFile,const string & _OrbFile,System_TzType * Ptr ):OneBodyFile(_OneBodyFile),OrbFile(_OrbFile),pSystem_Tz(Ptr),hbar_omega(Ptr->hbar_omega){}


  void setup(const string & _OneBodyFile,const string & _OrbFile,System_TzType * Ptr)
  {
    OneBodyFile=_OneBodyFile;
    OrbFile=_OrbFile;
    pSystem_Tz=Ptr;
    hbar_omega=Ptr->hbar_omega;
    setup();
  }
  
  /// setup orbitals,2B states and 2B mat.
  void setup()
  {
    setupOrbitals();
    setupTwoBodyStates();
    //setupTwoBodyMat();
  }
  void setupOrbitals();
  void setupTwoBodyMat();

  void get2BmatOnNeed(int a,int b,int c,int d,int J,TwoBodyMatEle & MatEle) const;

  void setupOverlap();
  
  void printOrbitals(int num=-1) const;

  complexd InOverOut(int orb) const;
  
  ///interface for HF cal.
  complexd get1B(int A,int bra,int ket) const;
  ///interface for HF cal.
  complexd get2B(int A,int a,int b,int c,int d,int J) const
  {
    TwoBodyMatEle MatEle;
    //    get2BmatAt(a,b,c,d,J,MatEle);
    get2BmatOnNeed(a,b,c,d,J,MatEle);
    if(MatEle.empty()) return 0;
    if(a==b) MatEle*=sqrt(2.0);
    if(c==d) MatEle*=sqrt(2.0);
    return MatEle[0]-pSystem_Tz->hbar_omega/A*MatEle[3];
  }

  ///interface for nocore shell model
  complexd get1B_p2(int bra,int ket) const;
  
  double hbar_omega;
  System_TzType *pSystem_Tz;

  string OrbFile;
  string OneBodyFile;
};
#endif
